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Abstract 


The  fundamental  problem  of  smoothing  and  differentiating  of  noisy  images  has  been 
previously  approached  in  two  different  ways:  1)  Minimization  of  a  smoothness  tunc- 
tional,  a  theoretically  well  understood  procedure  but  one  that  involves  the  solution  of 
a  very  large  system  of  equations  involving  all  the  pixels  of  the  image,  for  each  image. 
2)  Use  of  small  scale,  ready  made  filters  for  local  smoothing.  The  process  is  computa¬ 
tionally  cheap  but  generally  ad  hoc  and  not  very  reliable.  This  paper  offers  a  way  to 
combine  the  advantages  of  the  two  approaches.  We  construct  general  filters  for  local 
windows  of  the  image,  derived  from  maximization  of  smoothness  or  “regularization1' 
theory.  In  this  way  the  theoretically  robust  minimization  process  becomes  suitable  for 
practical  implementation,  possibly  in  real  time,  and  is  readily  adaptable  to  local  image 
properties.  Filters  for  more  reliable  derivatives  are  also  being  derived. 


Keywords:  smoothing,  surface  fitting,  differentiation,  interpolation,  minimal  curva-,  , 
ture,  filters,  splines,  regularization,  optimal  shape 
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1.  Introduction 


Smoothing  of  noisy  images  has  preoccupied  research  in  vision  and  image  processing  ever 
since  its  beginning,  and  a  variety  of  methods  have  been  developed  that  partly  overcome  the 
problem  but  leave  much  to  be  desired.  Strictly  speaking  the  problem  is  underdetermined 
since  there  is  no  way  to  distinguish  with  certainty  between  data  and  noise,  and  to  inter¬ 
polate  between  data  points,  so  some  reasonable  assumption  must  be  made.  Of  course  a 
theory  with  weak  extraneous  assumptions  is  preferable,  but  one  often  finds  it  unavoidable 
to  make  rather  strong  simplifying  assumptions  in  order  to  restrict  the  solution  space  and 
improve  computational  efficiency.  Thus  with  current  methods  one  has  to  trade  off  gener¬ 
ality  for  efficiency.  In  this  paper  we  show  that  one  can  relax  the  need  for  this  trade-off  to 
a  large  extent  and  achieve  a  reasonable  degree  of  both  generality  and  efficiency. 

We  first  brief  y  review  two  approaches  that  represent  different  choices  in  this  trade-off 
and  then  present  c  theory  that  combines  them. 

The  minimum  assumptions  approach  consists  of  a  minimization  of  some  cost  (or  en¬ 
ergy)  functional  E,  which  is  a  function  of  the  unknown  smoothed  image  /  and  the  data  g. 

( f,g  may  be  grey  levels,  depths,  or  any  interesting  function  defined  at  image  coordinates 
x.)  The  functional  can  be  expressed  generally  cis 

E  =  j  V(/,5)+  I  S(/) 

The  first  term  is  a  metric  V  that  measures  the  distance  between  the  desired  (minimizing) 
shape  /  and  the  data  g\  the  smaller  the  better.  The  second  term  integrates  a  smoothing 
function  S(f),  which  contains  derivatives  or  curvatures  of  the  shape  /.  This  generalizes 
the  idea  that  Grimson  [1981]  called  informally  “no  news  is  good  news”,  i.e.  the  surface’s 
departure  from  smoothness  is  minimal,  subject  to  the  requirement  of  small  distance  from 
the  data.  The  exact  forms  of  V',  S  depend  on  the  theory.  The  energy  minimization  principle 
was  used  by  Horn  [1983]  for  curves  and  by  Harrow  and  Tenenbaum  [1981]  for  interpreting 
line  drawings.  Poggio  et  al.  [1987]  have  put  this  principle  in  the  more  rigorous  mathe¬ 
matical  framework  of  the  “regularization  theory”  of  Tikhonov  and  Arsenin.  This  theory 
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shows  that  indeed,  under  weak  and  general  assumptions,  the  smoothing  term  f  S  turns 
an  ill-posed  problem  into  a  well-posed  one  in  the  sense  of  Hadamard  ( [  1923 j ) ,  i.e.  an 
underdetermined  problem  turns  into  one  that  has  a  unique  and  stable  solution.  Blake 
and  Zimmerman  [1987]  studied  the  convergence  and  other  properties  of  various  cost  func¬ 
tions,  Weiss  [1986]  investigated  cost  functions  that  adapt  themselves  to  local  properties  of 
the  shape,  and  Terzopoulos  [1988]  studied  implementation  methods  based  on  this  general 
framework. 

In  implementing  the  method,  each  pixel  is  assigned  an  unknown  variable  (or  several 
unknowns)  such  as  the  depth  of  the  image  at  that  point;  the  cost  function  is  calculated 
in  terms  of  the  unknowns.  Minimizing  the  cost  function  leads  to  equations  for  these 
unknowns.  In  this  way  one  obtains  a  large  set  of  coupled  equations  involving  all  the  pixels. 
Solving  the  system  for  3-D  shapes  has  shown  the  validity  of  the  basic  approach  but  also 
tne  prohibitive  cost,  even  with  acceleration  methods  such  as  multiple  grids.  Furthermore, 
since  the  solution  has  to  be  done  for  the  whole  image  at  once,  numerical  problems  arising 
from  a  small  area  can  impede  the  whole  solution.  Obviously  one  would  prefer  a  more 
local  solution,  that  involves  fewer  pixels  at  a  time  and  is  not  sensitive  to  every  numerical 
ill-conditioning  arising  from  anywhere  in  the  image. 

The  other  approach  to  smoothing  uses  filters  that  are  calculated  in  advance  for  a 
general  image,  and  are  used  by  convolving  them  with  the  particular  image  at  hand.  The 
filter  is  usually  only  a  few  pixels  wide,  making  the  convolution  efficient.  Gaussian  filters 
or  derivatives  of  them  (Marr,  [1989])  are  the  most  widely  utilized.  Among  many  other 
examples,  filters  utilizing  basis  functions  have  been  used  by  Hueckel  [1973]  and  Hummel 
[1979]  in  connection  with  edge  detection  (i.e.  detecting  non-smoothness),  and  by  Haralick 
[  1 984 ]  for  derivatives.  Recently,  Meer  and  Weiss  [1989]  derived  polynomial  filters,  based  on 
a  least  square  method,  for  smoothing  and  differentiation  to  high  order.  Another  group  of 
local  filters  is  based  on  statistical  “robust  estimation”  (Besl,  [1988]).  While  they  are  very 
good  at  rejecting  outliers,  they  are  not  always  as  good  in  dealing  with  the  rest  of  the  data, 
and  for  smoothing  images  with  relatively  little  hoise  they  are  sometimes  even  worse  than 
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the  Gaussian.  Besides,  these  filters  are  computationally  intensive,  requiring  an  iterative 
process. 

While  useful  in  many  applications,  many  of  the  local  filters  are  based  on  rather  ad  hoc 
assumptions,  so  it  is  not  surprising  that  their  output  is  sometimes  too  sensitive  to  both 
the  input  data  and  the  particular  details  of  the  filter  used.  Often  they  contain  internal 
parameters  that  need  to  be  carefully  fine  tuned  for  each  application  in  order  to  achieve 
a  result  that  looks  reasonable  to  humans.  It  is  clear  that  the  need  exists  for  local  filters 
based  on  a  solid  theoretical  foundation  such  as  the  regularization  theory  discussed  earlier. 

Calculating  derivatives  has  always  been  a  desirable  but  very  difficult  goal  in  vision. 
Global  regularization  theory  could  be  used  for  that  purpose,  but  then  the  computational 
cost  would  increase  even  more.  The  available  local  filters  have  proven  even  more  unreliable 
for  derivatives  than  for  smoothing.  Our  method  applies  differentiation  by  combining  the 
reliability  of  the  regularization  method  with  the  efficiency  of  the  local  filters. 

We  combine  the  two  methods  in  the  following  way:  we  define  a  minimal  curvature  cost 
function  and  apply  it  to  a  basic  “unit”  image,  i.e.  an  image  with  the  value  1  at  its  center 
and  0  elsewhere.  The  continuous  smoothing  functional  is  translated  to  the  discrete  mesh 
using  the  “finite  element”  method,  with  bivariate  Hermite  splines  as  the  basic  interpolation 
functions.  Minimizing  the  cost  function,  we  obtain  a  smoothed  version  of  the  unit  image. 
We  use  this  minimal  curvature  version  as  a  filter  or  mask  to  be  convolved  with  any  given 
image.  We  show  that  this  intuitive  filtering  process  can  be  justified  mathematically  by 
using  a  discrete  equivalent  of  the  Green  function  method  used  in  solving  differential  equa¬ 
tions.  While  most  of  the  treatment  involves  the  discrete  case,  we  do  study  some  properties 
of  the  continuous  analog. 

While  these  filters  have  merits  in  their  own  right,  they  can  be  of  even  greater  useful¬ 
ness  within  a  larger  context.  The  above  treatment,  like  other  treatments  throughout  the 
literature,  assumes  that  the  cost  function  has  constant  parameters.  In  a  common  physics 
analogy,  the  cost  function  is  likened  to  the  bending  energy  of  a  flexible  plate  that  represents 
the  surface.  These  methods  fails  when  a  surface  has  discontinuities,  or  sharp  edges.  A 
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usual  remedy  is  to  break  the  surface  into  smooth  parts  divided  by  the  discontinuities,  but 
this  immediately  raises  the  question  of  what  is  a  discontinuity  and  how  to  determine  it. 
Various  ad  hoc  and  not  very  satisfactory  attempts  have  been  made  to  solve  the  problem. 
Our  further  goal  is  to  solve  the  problem  by  varying  the  parameters  of  the  cost  function 
in  an  adaptive  way  that  will  fit  the  local  properties  of  the  surface.  A  surface  with  rapid 
changes  in  its  curvature  needs  a  more  flexible  “plate”  to  describe  it  than  does  a  slowly 
changing  part  of  the  image.  By  changing  the  physical  parameters  of  the  plate  according 
to  the  local  needs  of  the  surface  we  can  obtain  good  smoothing  properties  in  every  part  of 
the  surface  without  having  to  decide  prematurely  what  is  a  discontinuity. 

The  problem  was  hard  enough  to  solve  with  constant  parameters  because  of  the  large 
system  of  linear  equations  that  it  involved  Non-constar.t  parameters  wouid  a;  dee  the 
problem  nonlinear  and  almost  beyond  reach.  However,  our  reduction  of  the  solution  to  a 
simple  convolution  with  filters  not  only  makes  the  linear  solution  very  easy,  it  also  opens 
the  way  to  a  nonlinear  solution.  All  we  have  to  do  now  is  to  vary  the  parameters  of  the 
filter  during  the  convolution  in  accordance  with  the  local  properties  of  the  surface  we  deal 
with.  (In  a  way,  it  is  a  “context-dependent”  operation.)  In  a  subsequent  paper  we  will 
examine  various  ways  of  relating  the  filter  parameters,  or  the  flexibility  of  the  smoothing 
plate,  to  the  local  average  curvature  of  the  surface. 

In  the  next  sections  we  calculate  the  cost  function  as  a  function  of  image  variables 
by  using  1-D  basis  functions,  substitute  particular  spline  basis  functions,  derive  general 
solutions  and  put  them  in  the  form  of  a  convolution. 

2.  The  Cost  Function 

In  this  section  we  define  the  cost  function  to  be  minimized  and  express  it  in  terms  of 
unknown  image  variables.  This  is  done  by  expressing  the  2-D  image  by  general  basis 
function,  whose  coefficients  are  the  unknowns. 

The  function  to  minimize  is 

E  =  /  |  -  g[x))2dx  +  J  /n(V*/«  +  »** fly  +  »“f2n)dx  (1) 
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where  /  is  the  smoothed  (unknown)  image,  g  is  the  data  and  w  are  weights. 

The  first  term  is  simply  a  measure  of  the  distance  between  the  data  and  the  smoothed 
image,  which  is  a  quadratic  norm  with  weights.  The  second  part  is  the  smoothing,  or 
regularizing  term.  This  is  the  term  used  in  the  Tikhonov  and  Arsenin  theory  to  stabilize 
the  ill  posed  problem,  and  can  be  understood  by  a  physical  analogy  to  a  thin  plate.  The 
smoothing  term,  with  its  second  derivatives,  would  represent  a  bending  energy  term  of  the 
plate,  which  we  wish  to  minimize  along  with  the  distance  of  the  plate  from  the  data.  The 
coefficients  Aa,  with  a  =  xx,xy,yy,  are  the  regularization  coefficients  which  determine  the 
amount  of  smoothing.  The  bigger  the  smoothing  term,  the  smoother,  but  less  conforming 
to  the  data,  the  output  image  will  be.  If  A11  =  \yy  =  \xy /2  one  has  the  simple  isotropic 
bending  plate  (Landau  and  Lifshitz,  [1959]).  In  our  subsequent  work  we  will  be  interested 
in  AQ,u>  which  are  functions  of  both  the  coordinates  and  the  direction. 

We  now  have  to  discretize  the  problem.  We  could  represent  the  derivatives  in  the 
smoothing  term  by  finite  differences,  but  our  experience  (Weiss,  [1986])  and  others’  shows 
that  the  finite  element  method  is  more  stable,  besides  producing  a  smooth  shape  even 
between  the  pixels.  The  idea  is  simply  to  represent  the  surface  /(x)  as  a  linear  combination 
of  basis  functions  $,(x): 

m  =  £/.*,(*)  (2) 

t 

where  /,  are  unknown  coefficients  and  $,-(x)  are  given  basis  functions.  The  cost  function 
can  now  be  written  as 


E  =  J  j  YlltZ  Aa(x)(/,<?Q$,)2a!x  +  j  J  X)t£/(x)(/,-$,-(x)  -  g(x))2dx 

t  c*  i 

Minimizing  with  respect  to  a  particular  coefficient  we  obtain 

//?'■!?  Aa(dQ$t)(dQ$j)  +  |  dx  =  J  J  wg$jdx 


This  can  be  written  in  matrix  form  as 


KijU  =  Gj  (3) 
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where  the  vector  Gj  represents  the  image  data 

Gj  =  J  J  w(x)s(x)$,-(x)<foc 

and  Kij  is  the  analog  of  a  “stiffness  matrix”  of  the  plate: 


Kii  = 


(4) 


The  first  term  above  can  be  called  the  “mass”  term  K°° ,  and  the  others  are  the 
bending  terms,  Ka.  With  obvious  definitions  one  can  now  write 


Kij  =  K00  +  Kxx  +  Kxy  +  Kyy  (5) 

The  goal  of  finding  the  unknowns  /,  can  now  be  accomplished  by  solving  the  linear 
system  of  equations  (3)  for  the  unknowns  /,.  For  a  large  system  this  is  easier  said  than 
done  and  we  will  later  describe  an  efficient  direct  method  of  solution.  This  method  can  be 
applied  directly  in  the  case  of  sparse  data,  in  which  the  number  of  unknowns  is  relatively 
small. 

3.  Reduction  to  1-D 

So  far  we  ha,7e  dealt  with  general  basis  functions.  Since  all  functions  in  2-D  can  be 
expressed  by  basis  functions  that  are  products  of  1-D  functions,  we  will  express  KX]  in 
terms  of  such  simpler  functions.  Thus  we  introduce  the  1-D  functions  <f>k{x),4>k{y)-  For  the 
decompos^ion  one  has  to  use  different  i,j  indices  for  the  x  and  y  functions,  i.e.  iz ,  iy,Jz,Jy 
The  indices  of  the  stiffness  matrix  i,j  can  be  expressed  as  some  convenient  functions  of 
these  1-D  indices: 

i  —  *(**»  iy) 
j  =  j(jz, jy) 

fx  =  f\z  t  y  ’ 

$t(x)  =  4>lz{x)<t>ty(y), 


fj  —  fj.jy 
^(x)  = 


(7) 
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We  will  be  able  to  express  the  above  stiffness  matrix  terms  in  terms  of  three  basic  1-D 
matrices,  that  need  to  be  calculated  only  once  for  any  given  basis: 

K?t}x  =  J  4>ix(x)<f>Jz{x)dz 

Klz]z  =  J  dz(f>iAx)dx4>jI{x)dx  (8) 

Kfz]z  =  J  dxx4>iz{x)dxx4>ji{x)dx 


The  matrices  involving  y  are  similar. 

We  will  now  express  the  terms  X00,  Ka  with  a  =  xx,  xy,  yy  in  terms  of  the  above  basic 
matrices.  The  coefficients  w,  X  have  been  general  so  far,  but  to  perform  the  decomposition 
we  have  to  make  simplifying  assumptions  about  them.  We  assume  that  they  are  constant 
within  the  common  support  of  We  will  later  specialize  to  highly  localized  basis 

functions  of  the  size  of  a  pixel,  so  the  approximation  will  be  justified.  These  coefficients 
become  functions  of  i,j  rather  than  x  so  they  can  be  put  outside  the  integration. 

We  thus  have  for  the  mass  term 


Kif  =  /  J  w{x)4>iz  ix)<t>iv{y)<PjI  (x)<?>yy  {y)dxdy 

=  u>ij  J  <t>iz[x)<j>jz{x)dx  j  4>iy{y)<f>jv  [y)dy 


For  bendincr  terms  we  have,  first  for  the  x  direction 


K?/  =  J  J  ^xx (x)(3Zx<f)tI (x)4>iv [y)){dXx(pjz  (x)<f>jy  [y))dxdy 

=  *?*  J  dzx<t>tz{x)dxx<t>]z{x)dx  J  4>iy{y)(j>)y  ( y)dy 


(96) 


and  similarly  for  y  derivative 


(9c) 
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For  the  mixed  term 


K 


-  j  y  AIy(x)(5iy^,I(x)^,w(y))(ar!/^Ji(x)^s((y))dxdy 

=  Ai/  J  dT<{>tt{x)d^}Ax)dx  J  dv4> «„(y)<5y^j„(y)^y 


K} 


y  Jy 


(9d) 


Substituting  these  expressions  in  the  expression  for  the  total  stiffness  matrix,  eq.(o), 
completes  the  decomposition  into  1-D  quantities.  In  the  next  section  we  will  specialize  to 
specific  basis  functions. 


4.  The  Basis  Functions  and  the  Stiffness  Matrix 

The  general  expressions  for  the  stiffness  matrix  found  in  the  previous  section  can  be  used 
with  a  variety  of  basis  functions.  Meer  and  Weiss  [1989]  have  used  various  orthogonal 
polynomials  for  the  special  case  of  Aa  =  0,  i.e.  no  smoothing  term.  This  case  amounts  to  a 
simple  least  square  fitting,  and  it  was  easily  solved  analytically  because  the  stiffness  matrix 
is  diagonal.  We  obtained  good  smoothing  properties  but  the  derivatives  were  not  as  stable 
as  we  would  like.  Dierckx  [1977]  used  cubic  splines  for  the  above  special  case,  and  obtained 
a  system  of  equations  that  needs  to  be  solved.  In  this  paper  we  use  spline  functions  with 
narrow  local  supports,  i.e.  finite  elements,  to  make  the  solution  more  flexible  locally.  We 
then  eliminate  the  need  to  solve  a  system  of  equations  by  deriving  filters. 

Because  of  the  local  support  of  the  finite  elements,  the  geometry  of  the  image  points 
plays  an  important  role.  There  are  many  off  the  shelf  programs  for  solving  engineering 
problems  with  finite  elements.  In  a  general  finite  element  code  one  has  geometrical  “nodes" , 
i.e.  points  for  which  data  is  available  or  calculated.  They  are  input  annually  so  that  a 
wide  range  of  complicated  geometries  can  be  handled.  We  are  interested  here  in  a  simple 
geometry  that  can  be  dealt  with  automatically,  so  we  specialize  to  a  squait  domain  with 
a  rectangular  grid,  the  nodes  being  image  pixels.  Our  cost  function  is  somewhat  different 
from  what  can  be  found  in  the  ready  made  codes.  For  these  reasons,  after  experimenting 
with  some  of  these  programs  (ANSYS,  ADINA)  we  were  led  to  writing  our  own  program. 


The  stress  matrix  elements  found  from  (9)  all  involve  products  of  two  2-1)  basis 
functions,  each  of  which  consists  of  two  1-D  splines,  and  each  node  has  several  of  these 
functions,  so  a  good  bookkeeping  system  becomes  important. 

Since  the  K.  o  functions  are  nonzero  only  in  neighborhoods  of  pixels,  it  s  convenient  to 
enumerate  them  in  accordance  with  the  enumeration  of  the  pixels.  The  pixels  are  arranged 
on  _i.  rectangular  grid  and  are  denoted  by  mx,  my  (going  from  1  to  maxima  of  Mx ,  My.)  At 
each  node  (pixel)  we  have  a  small  number  of  basis  functions  so  the  indices  sz,sy  are 
used  to  distinguish  among  basis  functions  belonging  to  the  same  pixel.  These  two  indices 
enumerate  the  1-D  spline  functions  whose  products  make  up  the  2-D  basis  functions.  Each 
basis  function  has  two  equivalent  sets  of  indices:  While  the  i  (or  j)  indices  enumerate  all 
the  functions  in  a  linear  list  and  are  used  for  the  rows  and  columns  of  the  stiffness  matrix, 
the  m,s  indices  reflect  the  geometry  of  the  problem. 

A  basis  function  can  now  be  written  as  or  as  $>mi  ,my,St  ,«y  °r  for  brevity  as 
Similarly  for  its  coefficient  /,.  The  1-D  functions  can  be  written  as  <J>mz,sx,  4>my,ey- 

Fig.  1  shows  a  convenient  way  to  relate  to  two  kinds  of  indices  for  a  2  x  2  pixel 
example.  More  generally  we  can  establish: 

mx  (L  { 1 ,  Mz }  node  coordinate 

my  t  { 1 ,  A/y  }  node  coordinate 

m  -  (my  l)A/r  +  rnz  node  number 

,sx,.sy  r  {0,1}  intra  node  indices  (10) 

•s  2s T  -  sy 

i  -  4(m  1)  fs-M  -  4(my  -  l)Mz  +  4mx  +  2sz  +  sy  1 

j  ~  4(m,  1)  1  s'  f  1  -  4(my  -  l)Mx  +  4 m'x  *-  2sx  i  s'y  -f  1 

(We  note  that  m',$'  are  not  new  indices,  they  are  used  to  determine  the  column  j  of  the 
stiffness  matrix  while  m,. s  determine  i.) 
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One  can  get  a  good  grasp  of  the  basic  process  described  here  without  a  detailed 
understanding  of  the  plethora  of  indices  used,  but  the  intention  here  is  to  have  enough 
detail  for  an  actual  implementation. 

We  now  introduce  the  basis  functions  themselves.  Taking  the  distance  between  pixels 
to  be  1,  we  define  two  basic  functions  as  zero  outside  the  interval  [  —  1,1],  while  inside  they 
are  the  tv'o  cubic  Hermite  spline  polynomials  (Fig  2): 

<t>o,o(x)  =  (|x|  -  1)2(2|x|  +  1) 

,  (n) 

0O,l(x)  =  x(jxj  -  l)2 

To  form  the  whole  basis  these  are  shifted  so  they  are  centered  around  each  pixel,  in  both 
the  x  and  y  directions,  i.e. 

<£m,,ai(x)  =  <£o,8i(x  -  mx)  (12) 

/  —  fi{mzmysIB^)4’mz  a*  (y)  (13) 

t 

In  this  way  we  get  functions  with  local  support  in  a  rectangular  neighborhood  around 
each  pixel,  extending  only  up  to  the  next  pixel.  The  Hermite  polynomials  are  known  to 
have  good  properties  of  convergence  and  reliability  of  interpolation  up  to  second  order 
derivatives  (Strang  and  Fix,  [1973]). 

The  first  of  the  functions  in  eq.  (ll)  has  a  value  of  1  and  a  zero  derivative  at  the 
center  of  the  interval,  while  the  second  has  a  value  of  zero  and  a  derivative  1.  At  the  ends 
of  the  interval  both  the  values  and  the  derivatives  vanish,  so  there  is  no  (direct)  influence 
on  the  value  of  the  neighbors.  Thus  the  coefficients  /,  in  the  sum  (13)  are  in  fact  values 
of  the  function  /  ar.  ’.  its  derivatives.  It  is  easy  to  see  that 

/i(m,m,,0,0)  =  f{mx  i  my) 

/t  (mx  ,  1 ,0)  — 

(M) 

/i(m,m,l0,l)  —  dyf[mz,Tny) 

/i(mrmy,i,i)  dzdyf{mz ,  my) 

For  higher  derivatives,  higher  order  splines  are  needed  and  more  of  them  to  a  node,  but 
the  local  character  will  not  change. 
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Having  taken  care  of  the  bookkeeping  of  the  nodes  and  associated  basis  funct.ons,  we 
now  have  to  perform  the  integrations  of  eq.  (8).  Each  element  of  these  matrices  involve 
an  integration  over  a  pair  of  basis  functions  with  a  common  support.  Since  the  support 
is  local,  only  a  few  of  the  pairs  will  have  a  common  support.  Wc  can  divide  the  image 
intc  ectangular  tiles,  with  a  ncde  at  each  corner.  Each  tile  supports  only  a  small  number 
of  the  basis  functions.  Thus  one  .-an  save  effort  by  going  over  all  tiles  and  collecting  the 
contribution  of  each  to  a  matiix  element  of  Kn,  rather  than  going  over  all  pairs  i,j.  This 
is  known  as  the  “assembly”  process.  (Fix  and  Strang,  [1973],  Bathe,  [1982].)  From  the 
bookkeeping  point  of  view,  some  indices  are  replaced  by  ones  with  a  much  more  limited 
range: 

-  ruj  r  Am,,  mx Am'x 
my,m'y  — '  my  +  Amy,  my  -j-  A m'y 

where  mx,my  are  the  coordinates  of  the  bottom  left  node  and  Amx,  Amy  are  0  or  1 .  The 
general  indices  i,j  can  be  derived  from  them  with  eq.  (10): 

»  =  »(m*  +  Amz,  my  4-  Amy,sz,sy) 

3  -■  j{mx  +  A m'z,  my  +  A m'y,s'x,s'y)  (’5) 

To  assemble  the  matrices  we  now  go  over  the  indices  mz,  my,  A mx.  Amy,  A m'x,  A my,  s,  s' 
instead  of  mx,  my,  mx,  my,  s,  s'. 

Furthermore,  since  all  tiles  are  similar  it  is  enough  to  integrate  over  one  tile  and  use 
the  results  in  the  others.  The  actual  integration  can  be  done  in  1-E.  We  use  the  one-tile 
1-D  indices: 

ix  =  2A mx  +  sx  4-  1 
jx  —  2Amr  4  Oj  4  1 
iy  --  2A my  4-  sy  f  1 

jy  =  2Am'  J-  s'y  4  1  (16) 
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which  run  from  1  to  4.  The  factor  2  above  come  from  the  fact  that  there  are  currently  two 
basis  functions  per  node. 

Performing  the  integrations  (8)  on  the  one-tile  1-D  case  we  obtain  the  4x4  symmetric 
matrices,  for  an  interval  of  length  h 


Kn  =  I"  *.*3. 

Jo 


/ 156  22 

54 

-i3\ 

h 

4 

13 

-3 

420 

156 

-22 

=  f 


dxrf'ti 


30k 


dix<f>izdxx<f>]z  —  ^3 


4  J 

( 36  3  -36  3 

4  -3  -1 

36  -3 

4 

/  12  6  -12  6  \ 

4-6  2 

12  -6 

V  4  J 


and  similarly  for  y. 

The  above  quantities  are  the  contributions  from  each  tile  to  the  matrix  elements  of 
Kn  (n  =  1,2,3): 

=  Ku  (i- 

and  the  addition  to  the  total  stiffness  matrix  is  now  given  by  eqs.  (5,9)  as 

AATi,,  KlSxi<lu  +  1.K  l,+W„mMs.  kl b,+  Am.  ,m.K}.Ks, 

(18) 

These  terms  are  accumulated  in  the  general  stiffness  matrix  KX] .  There  is  usually  a 
contribution  from  more  than  one  tile  to  each  element  of  fft;. 

We  now  summarize  the  process  of  assembling  the  stiffness  matrix  as  follows: 
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For  all  image  pixels  mx,mv,  for  all  Amx,  Amy,  Am'x,  Am'y  equal  0  or  1, 
and  for  all  sx,  $y,  s'x,  equal  0  or  1  do: 

Calculate  the  one-tile  indices  ix,52,ty,iy,  eq.  (16). 

Retrieve  the  appropriate  element  of  the  4x4  matrices  if?  ^,...etc. 

Calculate  the  stiffness  matrix  indices  i,j,  eq.  (10). 

Calculate  the  contribution  AKXJ  from  eq.  (18)  and  add  to  the  appropriate 
array  element  Ktj. 

End  do 

This  stiffness  matrix  contains  all  our  assumptions  about  the  geometry  and  smoothness 
of  the  problem.  We  are  left  with  the  algebraic  problem  of  solving  the  system  of  equations 
that  it  represents,  eq.  (3).  This  is  dealt  with  next. 

5.  Reducing  the  Solution  to  Small  Filters 

In  theory,  solving  the  system  of  equations  (3)  should  be  a  straightforward  application  of 
algebraic  methods.  However,  an  image  with  1024  x  1024,  or  106  pixels  generates  a  stiffness 
matrix  with  1012  elements,  and  it  would  be  a  great  waste  of  resources  to  deal  with  it  in  a 
straightforward  way  even  if  it  could  be  done,  which  is  doubtful.  Fortunately,  the  matrix 
is  very  sparse,  since  only  a  few  of  the  pairs  of  basis  functions  have  a  common  support  (i.e. 
the  integrals  (8)  are  non-zero).  From  eq.  (10)  one  can  derive  the  pairs  t ,  j  that  have  a 
non-zero  element.  Previously  the  system  was  solved  using  iterative  methods  (Terzopoulos 
[1988]).  These  methods  are  very  efficient  in  use  of  memory,  but  many  iterations  over  the 
whole  image  are  needed  to  get  close  to  convergence.  Direct  methods  give  a  closed  form 
solution  but  may  need  sophisticated  hashing  techniques  to  save  memory,  depending  on  the 
structure  of  the  matrix.  (Ortega,  [1981].)  In  our  case  this  structure  is  quite  simple,  being 
banded  along  the  diagonal,  so  storage  is  quite  simple.  Direct  methods  are  now  used  in 
most  finite  element  codes  (Bathe,  [1982])  and  we  adopt  one  here. 

In  this  section  we  show  that  the  one  can  solve  the  problem  by  solving  for  small  windows 
of  the  image,  provided  that  we  find  a  general  solution  for  any  small  image.  For  this  process 
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a  direct  method  becomes  even  more  attractive,  as  the  small  window  makes  it  very  easy  to 
apply. 

We  begin  by  describing  briefly  the  continuous  analog  of  the  method.  We  are  given 
the  differential  equation,  which  we  keep  in  1-D  for  simplicity: 

Dxf{x)  =  g{x) 

where  Dx  is  a  differential  operator,  f(x)  is  an  unknown  function  and  g( x)  is  a  given  term, 
e.g.  data  (the  inhomogeneous  term).  Now,  instead  of  solving  for  each  g  one  can  solve  for 
particularly  simple  data,  i.e.  delta  functions,  and  build  a  general  solution  on  this  basis. 
Thus  we  define  the  “Green  function”  h(x,t )  as  the  solution  of 

Drh(x,t)  =  <5(z  -  t) 


The  general  solution  is  now  given  by 

f  ~  J  h(x,t)f(t)dx 


Proof: 

Dxf  —  j  Dxh(x,t)f{t)  —  j  6(x-t)f(t)=f(x) 

If  all  the  coefficients  in  the  differential  operator  are  constant,  then  the  equation  is  shift 
invariant,  i.e.  we  have  the  same  Green  function  around  each  point  t: 


h(x ,  t)  =  h(x  —  t) 


and  thus  the  general  solution  above  becomes  a  convolution. 

We  will  later  find  the  Green  function  for  the  continuous  version  of  our  smoothing 
problem. 

We  return  now  to  our  discrete  case.  We  have  (eq.  (3)): 

J 


14 


The  indices  i,j  are  functions  of  the  geometrical  node  numbers  m,s  (eq.  (10)).  We 


m  =  (mx,my)  s  =  s(sx,sy) 
i  =  i(m,s)  j-j(m',s') 

and  ^  enumerates  the  basis  functions  at  each  node.  So  we  can  write  eq.  (3)  as 

Km,m'  ,e,s'  fm',e'  —  Cm, 3 

m' ,3' 

We  can  simplify  matters  if  we  assume  that  the  data  is  concentrated  in  the  node,  a 
standard  practice.  From  eqs.  (4,7,11)  one  can  see  that  Gmig  vanishes  for  all  s  except  s  =  0, 
in  which  case  it  is  equal  to  the  given  data  at  the  pixel  m,  Gm  =  Gmio-  Thus  the  solution 
for  the  above  equation  vanishes  for  s  ^  0  and  we  can  limit  ourselves  to  the  solution  of 

^  '  Krn,m'  ,0,8*  fm'  ,s'  —  G  m  (19) 

m'  ,3' 

In  analogy  with  the  continuous  case  we  will  now  look  for  a  solution  to  a  system  of 
equations  with  the  data  being  a  Kronecker  delta  rather  than  a  general  G: 

^  ^  ^m,m' ,0,j'^(^  1  *  5  ft)  —  bmn  (20) 

m'  ,3' 

where  n  —  n(nx,ny )  is  an  arbitrary  node  of  the  grid. 

The  solution  of  (19)  for  general  data  is  now 

fm,a  =  ^2 


Proof: 


y  y  Km,m' , 0,s'  y  ^  h[m  ,S  ,  ft)  Gn  —  y  ^  fimnGn  —  Gr 


In  dealing  with  images  we  observe  that  our  problem  is  shift  invariant  in  the  node 
number  m,  because  we  constructed  the  finite  element  formalism  in  this  way,  i.e.  the 
same  basis  functions  just  shifted  unchanged  from  one  pixel  to  another.  (We  assume  for 
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the  moment  an  unbounded  image,  analogous  to  the  continuous  case.)  Thus  our  Green 
function  h  loses  one  variable  and  becomes  a  filter  h(m  —  n,s)  which  is  convolved  with  the 
image: 

fm,s  =  -  n,s)Gn  (21) 

n 

Furthermore,  it  is  reasonable  to  assume  that  the  quantities  at  pixel  m  will  not  be 
influenced  very  much  by  the  data  in  a  pixel  n  which  is  far  away  from  m.  Therefore  the 
summation  above  does  not  need  to  be  carried  out  over  the  whole  image  but  can  be  limited 
to  some  rectangular  window  extending  Nx,Ny  pixels  around  m.  Thus  we  can  finally  write 
the  solution  in  term  of  the  filters  as: 

fmxmvaxsv  =  ^  h{mx  —  nx,my  —  ny,  sx,  Sy)Gnjtrly  (22) 

nx=-Nx,ffz 

ny  —  —Ny,Ny 

Thus,  the  solution  at  the  pixel  mx,my  is  given  by  a  known  mask  applied  to  a  window 
around  that  pixel. 

One  can  see  from  (22)  that  there  is  a  separate  filter  for  each  sx,sy,  with  eq.  (14) 
giving  the  meaning  of  each  one.  Thus, 

h(mx,my,  0,0) 

will  yield  the  smoothed  function  /  itself,  while 

h(mx,my,  1,0) 

is  the  x  derivative,  and  the  combinations  (0,1),  (1,1)  yield  the  y  and  x,  y  derivatives 
respectively.  Since  convolution  is  associative,  we  can  concatenate  two  first  derivative  filters 
to  obtain  one  for  the  second  derivative: 

dXzf{rn)  ~  h(m,  1,0)  0  ( h(m ,  1,0)  <g>  /(m))  =  ( h(m ,  1,0)  0  ( h(m ,  1,0))  0  f(m) 

so  that  the  filter  for  the  second  derivative  is 

h(m,2,0)  =  h(m,  1,0)  0  h(m,  1,0) 
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and  similarly  for  other  derivatives.  There  is  no  point,  however,  in  continuing  the  process 
beyond  the  order  of  the  spline  used:  we  need  higher  order  splines  than  our  cubics  to  obtain 
third  and  fourth  derivatives. 

In  conclusion,  we  have  reduced  the  regularization  problem,  which  led  to  a  large  system 
of  equations,  to  convolutions  of  the  image  with  small  filters.  The  filters  can  be  derived 
once  and  for  all  (for  a  cost  function  with  given  factors)  by  solving  eq.  (20)  for  a  small 
window,  a  relatively  easy  task.  The  results  can  be  stored  for  use  with  any  input  images. 
The  next  section  summarizes  the  solution  method  used  for  the  window. 

6.  Solving  for  the  Filters 

Deriving  the  filters  amounts  to  solving  the  system  of  linear  equations  (20)  for  a  small 
window.  In  principle,  this  could  be  done  by  the  standard  method  of  Gaussian  elimination, 
but  then  one  has  to  deal  with  the  recurring  problems  of  ill-conditioning  and  numerical 
instabilities.  A  method  applicable  to  symmetric  positive  definite  matrices  is  the  Cholesky 
decomposition  (Ortega,  [1981]).  It  guarantees  well  conditioning  and  numerical  stability. 
We  describe  here  a  modification  for  a  general  symmetric  (non-singular)  matrix.  A  stiffness 
matrix  of  a  stable  physical  system  is  always  positive  definite.  In  our  case  it  is  at  least 
symmetric. 

We  start  with  the  decomposition  of  the  symmetric  matrix  into  a  product: 

K  =  LDLt 

where  L  is  a  lower  triangular  matrix  having  l’s  on  its  diagonal 


Lt  is  its  transpose,  an  upper  triangular  matrix,  and  D  is  a  diagonal  matrix.  Explicitly: 

n 

Kij  —  ^  LikDkLjk 

k= l 
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where  n  in  the  matrix  size,  in  our  case  the  total  number  of  unknowns  in  the  image.  We 
further  assume,  without  loss  of  generality,  that  the  matrix  is  banded  with  p  non-zero 
diagonals  above  and  below  the  main  diagonal.  This  is  the  case  for  our  sparse  stiffness 
matrix.  The  following  iterative  algorithm,  based  on  the  above  summation,  can  be  verified 
by  induction: 


For  j  =  1, . . . ,  n  do: 

q  ■>  -  iucix(l,j  p) 
i- 1 

”k  -  Ksj  -  £  L%Dk 

k=q 

For  i  ~  j  +  1,  •  •  • ,  min(y  +  p,  n)  do: 
r  <—  max(l,  t  —  p) 

y-i 

Lij  (Kij  -  £  KikKjkDk)/D j 

k=zr 

End  do 

End  do 

The  output  elements  L,y  also  form  a  banded  matrix  and  can  be  conveniently  stored 
by  overwriting  the  corresponding  input  elements  Kij  as  soon  they  are  calculated.  The 
diagonal  matrix  elements  Dk  can  be  stored  in  place  on  the  l’s  in  L.  They  also  happen  to 
be  the  eigenvalues  of  the  matrix. 

We  can  see  from  the  limits  of  the  iterations  and  the  summation  that  the  decomposition 
is  of  linear  complexity  in  n  and  0(p2)  in  the  bandwidth,  i.e.  it  is  0(np2).  For  a  non-banded 
matrix  we  have  q  =  r  =  1  and  the  complexity  would  be  n3.  Thus  if  the  bandwidth  is  much 
smaller  than  the  matrix  dimension  one  obtains  great  savings. 

For  given  data  G  the  system  of  equations 


Kf  -  G 


can  now  be  solved  in  two  stages,  with  an  intermediate  solution  vector  /',  as  follows  (“back 
substitution”): 

Lf  -  G 


DLTf  =  /' 


(24) 
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The  first  term  can  be  written  explicitly  as 

£w;  =  c. 

k 

As  before,  by  separating  out  the  last  term  of  the  sum  one  obtains  the  iterative  algorithm: 

For  i  =  1, . . .  ,n  do: 

r  max(l,  t  —  p) 
i  —  1 

/;  <-  <?i  -  E L 

h  =  r 

End  do 

Similarly,  the  second  stage  is 

DtJ2Lkifk  =  fl 

k 

yielding  the  “mirror”  algorithm 

For  t  =  n, . . . ,  1  do: 

r  <—  min(l,t  +  p) 

r  +  1 

fi-  f'jDi-Y,Lkifk 

k  =  r 

End  do 

For  a  banded  matrix,  the  two  stages  are  O(np)  in  complexity,  while  for  the  non- 
banded  case  (r  =  1)  we  have  0(n2).  Thus  the  main  effort  is  in  the  decomposition  stage, 
with  0(np2).  In  our  implementation  the  band  width  is  proportional  to  the  window  side 
N  and  the  matrix  dimension  to  N 2  so  we  have  a  total  of  0(N 4).  A  window  with  N  «  10 
takes  a  few  seconds  for  a  VAX  to  solve. 

Finally,  we  have  to  substitute  the  appropriate  data  G  in  the  above  algorithms.  Of 
course  the  algorithms  could  be  used  with  the  whole  image  itself  as  data.  That  would  be 
useful  for  sparse  data  with  a  limited  number  of  unknowns.  In  deriving  the  filters,  however, 
the  data  is  “unit”  images  (Kronecker  <5s)  according  to  eq.  (20).  Denoting  the  central  pixel 
in  a  window  by  nz,ny  we  want  to  solve  the  equation 

fi(m,,mv,a)  =  ^mtn,^myny 
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The  solution  can  be  broken  into  S  parts,  one  for  each  s  =  1, . . . ,  S,  corresponding  to  our 
5  filters  (22).  Algebraically,  we  have  obtained  one  row  of  the  inverse  of  the  matrix  and 
there  is  no  need  to  invert  the  whole  matrix. 


7.  The  Continuous  Approximation 

For  the  continuous  case  the  functional  (l)  can  be  minimized  using  the  Euler-Lagrange 
equations  of  the  calculus  of  variation,  yielding  the  fourth  order  differential  equation 

^V4$(x)+$(x)  =  <7(x) 

In  cylindrical  coordinates  the  rotationally  symmetric  Green  function  involves  JV0(\/fr),  i.e. 
Bessel  functions  of  the  second  kind  with  complex  variable.  The  function  diverges  at  zero 
and  its  asymptotic  behavior  at  infinity  is  unclear.  In  Cartesian  coordinates  we  can  learn 
the  main  feature  of  the  filters  by  examining  the  solution  of  the  1-D  equation 


4 

dx4 


+  < p(x )  —  6(x  —  t) 


where  we  have  defined  4a4  =  X/w. 
The  Green  function  is 


h{x  —  t) 


x  —  t 

cos - 

2a 


Substituting  h(x  —  t )  in  the  above  equation,  it  is  easy  to  verify  that  the  left  hand  side 
vanishes  for  x  ^  t,  and  that  its  integral  has  a  jump  of  1  over  the  discontinuity  at  x  =  t, 
(due  to  the  jump  in  the  third  derivative),  making  it  a  delta  function.  The  parameter  o 
was  chosen  so  that  the  Green  function  will  behave  generally  like  a  Gaussian  so  that  the 
two  can  be  compared.  Fig.  3  shows  the  Green  function  in  a  solid  line  and  the  Gaussian. 
Fig.  4  shows  their  respective  derivatives.  These  are  used  in  convolution  with  an  image 
to  obtain  an  image  derivative.  We  see  that  one  major  difference  is  in  the  central  parts 
of  the  filters.  The  bending  Green  function  is  less  smooth  at  the  top,  and  its  derivative  is 
discontinuous.  This  can  make  a  difference  for  features  smaller  than  a  or  high  frequencies. 
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Indeed,  the  Fourier  transform  of  the  filters  (Figs.  5,6)  shows  that  the  Gaussian  filters, 
particularly  the  differentiation  filter,  attenuate  higher  frequencies  more  strongly  than  the 
bending  filter.  Another  difference  is  the  “bump”  in  the  Fourier  transform  of  the  Green 
function,  resulting  from  the  cos(x/2a)  factor.  Thus  the  function  looks  somewhat  like  a 
band-pass  filter.  A  combination  of  such  filters  with  different  scales  a  may  be  tailored  to 
provide  a  flatter  band-pass  filter  between  given  frequencies,  something  that  can’t  be  done 
with  the  Gaussian.  (The  expressions  depicted  are  e-cr2/2  and  (2a2  +  1) / (4  a4  +  l).) 

8.  Experiments 

We  solved  the  linear  equations  (20)  for  various  window  sizes  and  convolved  the  resulting 
filter  with  real  images.  Fig.  7  shows  a  cross  section  through  the  center  of  the  smoothing 
filter,  and  Fig.  8  does  the  same  with  the  filter  for  the  first  derivative.  Fig.  9  shows  a 
smoothed  image,  and  Figs.  10,11,12  are  derivatives  dx,dy,dyy  respectively.  While  any 
effects  are  hard  to  measure  quantitatively  for  real  images,  we  can  see  that  the  derivatives 
look  much  like  one  would  expect.  The  horizontal  edges  almost  disappear  with  the  dx  filter 
and  they  show  progressively  more  variation  for  increasing  differentiation  order  dy,dyy. 
The  smoother  parts  are  progressively  flattened. 

Fig.  13  shows  an  isolated  edge  smoothed  by  the  Green  function  (solid  line)  and  a 
Gaussian.  The  Green  smoothing  correctly  represents  the  top  of  the  edge  but  needs  some 
distance  until  it  clings  to  the  straight  sides.  This  can  be  expected  from  a  flexible  plate 
that  one  tries  to  fit  over  the  edge.  The  Gaussian  smooths  the  edge  away,  consistent  with 
a  dissipating  temperature  gradient.  Clearly  both  filters  in  their  simple  form  would  better 
be  used  on  smoother  parts  of  the  image.  A  subsequent  paper  will  deal  with  cost  functions 
whose  factors  adapt  locally  to  the  properties  of  the  image,  so  that  a  will  get  smaller  as  the 
filter  gets  closer  to  the  edge. 

9.  Conclusion 

We  have  derived  filters  for  smoothing  and  differentiation  based  on  the  principle  of  minimal 
curvature  of  surfaces,  on  rectangular  grids.  This  principle  is  based  on  the  regularization 
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theory  that  provides  reliable  stabilization  of  noisy  or  insufficient  data.  The  derivation  was 
numerical  for  the  discrete  case  and  analytic  for  the  continuous  approximation. 

One  can  use  certain  parameters  to  control  the  degree  of  smoothing,  unlike  least  square 
fitting,  which  is  obtained  as  a  special  case.  The  filters  turn  out  to  have  some  general 
similarities  to  the  Gaussian  but  they  preserve  some  frequencies  better  than  the  Gaussian. 

The  filtered  data  can  be  used  both  as  image  values  and  derivatives  at  the  pixels, 
and  as  factors  that  multiply  spline  functions  interpolating  between  the  pixels.  The  latter 
capability  is  particularly  useful  for  sparse  data.  Other  methods  need  to  solve  a  large  system 
of  equations  for  these  interpolation  factors. 

For  irregular  data  with  varying  weights  no  simple  filter  method  is  applicable  and  a 
whole  system  of  linear  equations  must  be  solved.  The  direct  solution  method  we  described 
is  particularly  stable  and  efficient  for  this  purpose. 
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Figure  Captions 


Fig.  1.  Enumeration  schemes  for  local  basis  functions 
Fig.  2.  Cubic  Hermite  splines 

Fig.  3.  Smoothing  filters:  minimal-curvature  (solid  line)  and  Gaussian 

Fig.  4.  Differentiation  filters:  minimal-curvature  (solid  line)  and  Gaussian 

Fig.  5.  Fourier  transform  of  the  smoothing  filters 

Fig.  6.  Fourier  transform  of  the  differentiation  filters 

Fig.  7.  Discrete  minimal-curvature  filter  for  smoothing 

Fig.  8.  Discrete  minimal-curvature  filter  for  differentiation 

Fig.  9.  Smoothed  image 

Fig.  10.  x-derivative  of  image 

Fig.  11.  y-derivative  of  image 

Fig.  12.  Second  y-derivativc  of  image 

Fig.  13.  Edge  smoothing 
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